Lattice Boltzmann simulations of liquid crystalline fluids: 

active gels and blue phases 
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Lattice Boltzmann simulations have become a 
method of choice to solve the hydrodynamic 
equations of motion of a number of complex flu- 
ids. Here we review some recent applications of 
lattice Boltzmann to study the hydrodynamics 
of liquid crystalline materials. In particular, we 
focus on the study of (a) the exotic blue phases 
of cholesteric liquid crystals, and (b) active gels 
- a model system for actin plus myosin solutions 
or bacterial suspensions. In both cases lattice 
Boltzmann studies have proved useful to pro- 
vide new insights into these complex materials. 

1 Introduction 

In recent years, the lattice Boltzmann (LB) al- 
gorithm [TJ has emerged powerful method 
to study fluid dynamics. Due to its conceptual 
simplicity and codability (particularly on paral- 
lel computers), LB provides an attractive alter- 
native to other methods such as finite elements 
algorithms. In the last couple of decades, in 
particular, the LB method has increasingly been 
applied to the hydrodynamics of complex fluids, 
such as binary fluids, colloidal suspensions and 
liquid crystals [21 El HI El El IB] ■ 
In applying LB to complex fluids, one often aims 
at solving two coupled sets of partial differential 
equations. One set describes the evolution of 
the order parameter (e.g., composition for bi- 
nary mixtures, or orientational order in liquid 
crystals), whereas the other set describes con- 
servation of mass and momentum (via the con- 



tinuity and Navier-Stokes equations for the ve- 
locity field). In liquid crystals, which are the 
focus of this article, one typically considers the 
Beris-Edwards model [9] , which is defined start- 
ing from a free energy expressed in terms of 
a tensorial order parameter, Q, whose largest 
eigenvalue describes locally the strength of lo- 
cal molecular alignment (nematic order), and 
whose corresponding eigenvector defines the di- 
rector n along which this alignment prevails. 
(See the Appendix for the precise form of the 
free energy adopted.) Note that simpler descrip- 
tions, in which the magnitude of the ordering is 
assumed constant and only the director varies, 
are unsuitable for describing blue phases. This 
is because such phases contain defect lines at 
which the ordering drops locally [TO] . 
The equation of motion for Q is then [9]: 



D t Q = TH. 



(1) 



The left member of equation [I] is a "material 
derivative" describing the time evolution of the 
order parameter advected with the velocity u 
of a fluid element. For fluids of rod- like parti- 
cles such as liquid crystals, flow gradients may 
lead to local rotations and these are also taken 
into account by the material derivative. In Eq. 
(H r is a collective rotational diffusion constant, 
which sets the time scale for the relaxation of 
orientational order (usually in the millisecond 
range for small- molecule liquid crystals), and H 
is the "molecular field" , which provides the force 
for this relaxation motion. The molecular field 



1 



involves the derivative of the free energy with re- 
spect to the order parameter (see the Appendix 
for the specific form of H). 
As stated previously, the fluid velocity obeys the 
continuity equation, d t p = — V.u = 8 a u a where 
p is the fluid density, which in practical cases 
this can be taken as constant (so that the fluid 
is incompressible). However, a feature of LB 
is that slight fluid compressibility is maintained 
in the algorithm; this makes all the dynamics 
fully local (the sound-speed is finite) rather than 
having to solve each timestep for a pressure field 
that responds instantaneously to distant events. 
(Locality is important to efficient parallelization 
strategies, and maintains near-linear scaling of 
the computational cost with the size of the sys- 
tem under investigation [IT].) 
The fluid velocity also obeys the Navier-Stokes 
equation, which for effectively incompressible 
fluids reads 

p[d t + updp]u a = dpIL a p (2) 

+ 7)8/3 {daU/i + d/3U a ) , 

where r/ is the viscosity, Cartesian components 
are denoted by Greek indices, and is a ther- 
modynamic stress tensor. This tensor, like H, 
is found by differentiation of the free energy, 
and its divergence dpll a p represents an effective 
body force, acting on the fluid, arising from the 
response to deformation of the order parameter 
field. (Appendix A gives the explicit form.) In 
the case of active fluids, the forces do not solely 
stem from a free energy but include terms aris- 
ing from the dissipative conversion of chemical 
energy into motion. This creates an additional 
contribution to H a p. 

Eqs. [1] and [2] are very complex to solve due to 
their inherent nonlinearities, and only limited 
progress is possible with analytical techniques. 
In contrast, LB offers an ideal method to solve 
numerically these equations, allowing their full 
dynamics to be addressed not only in one di- 
mension but also in 2D and 3D flows. Such 
studies can not only test approximate analytic 
solutions where these exist, but also lead to im- 
proved insight into the nonlinear physics con- 
tained in the underlying models. 
For purely Newtonian fluids, the LB algorithm 
proceeds by the introduction of "mesoscopic" 



velocity distribution functions /(cjjx), propor- 
tional to the density of (notional) fluid particles 
sitting at a lattice node at x having a certain ve- 
locity Cj chosen from a discrete set pp. (These 
discrete velocities correspond to moving by one 
lattice site in a timestep.) The distribution 
functions evolve according to an appropriate lo- 
cal dynamics, which recovers the Navier-Stokes 
and continuity equations upon coarse graining. 
(The density p then equates to J2i f( c i) an d the 
average fluid velocity u becomes the first mo- 
ment, J2if( c i) c i, of the distribution.) 
Historically, the first LB approach to study com- 
plex fluids [31 [5] consisted of one extra set of dis- 
tribution functions for each new order parame- 
ter component entering the equations of motion. 
More recently, however, it has become apparent 
that this "full LB" approach had the drawback 
of requiring a large memory to store all the dis- 
tribution functions for the whole lattice, on top 
of being rather cumbersome from a theoretical 
viewpoint. At the same time, the new genera- 
tion of LB studies for complex fluids has shown 
the need to address ever larger systems, to make 
best use of the potentially very good scalability 
of parallel LB codes. 

As a result, new hybrid algorithms have been 
coded and deployed, for both binary fluids and 
liquid crystals [Hl[T3], where the LB algorithm 
is used to solve the Navier-Stokes equation (Eq. 
[2]), and is coupled to a standard finite-difference 
solver for the order parameter dynamics (Eq. 
[TJ. At each timestep, the fluid velocity found 
by LB is used to calculate the advective deriva- 
tive in Eq. [H while the order parameter found 
from that equation is used to compute the forc- 
ing term, dpIL a p, in Eq. [2J By this division of 
labour, LB is only used to handle the momen- 
tum and mass transport - the problem for which 
it was originally devised. 

In what follows, we review recent applications of 
full and hybrid LB algorithms to the study of ac- 
tive liquid crystalline fluids and the blue phases 
of cholesteric liquid crystals. In both cases, as 
we shall see, numerical simulations have proved 
extremely helpful in providing a link between 
theoretical predictions and experimental obser- 
vations. 

While in this work we focus on a "hydrody- 



2 



namic" description of liquid crystalline fluids 
via continuum models, it is important to note 
that there is a variety of other coarse grained 
models and methods to study liquid crystals. 
Most relevant to our topics here, liquid crys- 
talline molecules may be individually modelled 
via e.g. soft spherocylinders, interacting with 
Gay-Berne potentials (see e.g. [HI [15] for a se- 
ries of recent examples) and a related approach 
has also managed to stabilise blue phases [16] . 
although the length scale accessible with this 
more detailed approach is typically significantly 
smaller than the ones which can be studied with 
an LB continuum description. 

2 Active fluid simulations 

Active fluids have become a highly topical re- 
search area at the interface of soft matter and 
biological physics. Generally speaking, an ac- 
tive particle "absorbs energy from its surround- 
ings or from an internal fuel tank and dissipates 
it in the process of carrying out internal move- 
ments". (The quote is from a well-known pa- 
per of the Bangalore group, which helped pi- 
oneer the use of continuum models of active 
matter [IT].) This definition applies to bacte- 
ria, swimming algae and other microorganisms, 
as well as to living cells or cell extracts. How- 
ever active materials may also be non-biological, 
and a synthetic example is a shaken granular 
fluid [TB] . As the definition clearly implies, ac- 
tive materials may remain far from thermody- 
namic equilibrium even in steady state, due to 
their continuous energy intake, and this renders 
their properties of particular interest to physi- 
cists. 

We consider active fluids comprising a con- 
centrated suspension of such active particles. 
Paradigmatic examples of active fluids are bac- 
terial suspensions and solutions of cytoskeletal 
gel components (actin fibers or microtubules) 
with molecular motors (myosin or kinesin). The 
term "active gels" is also widely used for these 
materials; but while all are non-Newtonian, not 
all of them are strongly viscoelastic. Experi- 
mental studies with active fluids have uncov- 
ered a wide range of intriguing and non-trivial 



physical properties. For instance, microscopy 
studies of droplets of Salmonella [19] and of B. 
subtilis [20j EH [22] reveal that when the bacte- 
ria are concentrated enough (more than about 
20-30% in volume fraction), long range corre- 
lations arise, creating eye-catching patterns of 
flow involving long-lived vortices (see Fig. 1). 
These resemble the turbulent flow of fluids at 
high Reynold numbers, although remarkably in 
this case, the Reynolds number is very small - 
effectively zero [23] • (In this respect, the result- 
ing "bacterial turbulence" resembles elastic tur- 
bulence in polymer solutions, which is attained 
past a critical value of the Wiessenberg num- 
ber [23 ES].) 

The equations of motion of an active liquid crys- 
talline fluid have been written down, either on 
the basis of symmetry [T7J |25], or via a coarse 
graining of an underlying microscopic model of 
stiff cytoskeletal gels and molecular motors [2Z] ■ 
There are three main differences between these 
equations and the equations of motion of a pas- 
sive liquid crystal. The first one is the presence 
of an active term in the stress tensor, whose 
divergence acts CIS db force in the Navier- Stokes 
equation (Eq. [2]). This extra active term has 
been first shown in [2H] to be proportional to 
an activity constant, (, which to lower order is 
linear in the energy uptake of the fluid (e.g., 
via ATP hydrolysis), times the local order pa- 
rameter Q. A second term arises from activ- 
ity also, but the form of this can be absorbed 
into the free energy (see Appendix) and we set 
this to zero. A third and final difference can 
arise in cases where an oriented swarm of swim- 
ming particles have a collective mean velocity 
relative to the surrounding fluid, causing a 'self- 
advection' effect. Such materials are called 'po- 
lar' and are distinct from the 'apolar' case which 
describes either a swarm in which equal number 
of particles swim forward and backward along 
the director axis, or particles which are non- 
motile but nonetheless exert active forces on the 
fluid. (Motile and non-motile active particles 
are sometimes called "movers" and "shakers" 
respectively.) More details on the exact forms 
of these active terms, and of the equations of 
motion used to describe polar systems, are given 
in the Appendix. 
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The sign of £ - hence of the active apolar con- 
tribution to the stress tensor - is of vital im- 
portance for the hydrodynamics and the rhe- 
ology of active fluids [7J [T7]. A positive ( 
corresponds to a suspension of extensile active 
particle, or "pushers" [22], which exert forces 
along the molecular axis away from the cen- 
tre of mass and towards the surrounding fluid. 
A negative ( corresponds to a fluid of contrac- 
tile active particles, or "pullers" [29], for which 
the force dipoles is exerted axially towards the 
centre of mass. Examples of contractile fluids 
are suspensions of the biflagellated alga Chlamy- 
domonas and actomyosin gels (or more general 
suspensions of non-permanently cross-linked cy- 
toskeletal gels and molecular motors). On the 
other hand, the majority of bacteria are thought 
to be extensile [29]. However, we are currently 
lacking quantitative experiments which measure 
the velocity field around active particles, which 
could lead to estimates of the values of ( in these 
various cases. 

Alongside the two distinctions already made 
(apolar versus polar, and contractile versus ex- 
tensile), all rodlike molecules, including passive 
ones, fall into two further categories, known as 
"flow aligning" and "flow tumbling". With- 
out activity, the former exhibit stable flow in 
which the director is inclined to the flow direc- 
tion at a certain angle (the Leslie angle) whereas 
the latter undergo continuous director evolu- 
tion, which is frequently chaotic [20] • To avoid 
the double complexity created by the flow tum- 
bling instability on top of activity, we address 
here only the flow- aligning case. 
Early theoretical work determined the linear 
stability of these systems and it was found 
that an infinite sample of active material with 
nonzero Q (whether polar or apolar) is hy- 
drodynamically unstable to order parameter 
fluctuations, and that this instability is con- 
nected to the generation of a spontaneous fluid 
flow [13 ESI M\- (For extensile - but not con- 
tractile - particles, this instability is present 
even if the flow velocity is constrained to be 
a function of one spatial coordinate only.) It 
was also realised that the introduction of bound- 
aries together with suitable "anchoring condi- 
tions" (fixing the molecular orientation) would 



lead to the stabilisation of the non-flowing or- 
dered phase of uniform Q. (Stability is restored 
for small values of lying below a threshold ( c 
which decreases with system size, L, as \/L 2 .) 
However, within these analytical approaches, it 
was not possible to determine the ultimate flow 
pattern resulting from this hydrodynamic insta- 
bility. 

The simulations reported in [6l[7J[T3] gave there- 
fore the first quantitative predictions for the 
resulting spontaneous flow patterns in unsta- 
ble active fluids. We focus first on a quasi- ID 
slab geometry with planar anchoring along the 
boundary. Here it was found that for an apo- 
lar fluid, upon increasing the activity, the sys- 
tem organises into a spontaneous Poiseuille flow 
(with a smoothly varying flow velocity, maximal 
at the centre of the slab). This spontaneously 
breaks symmetry and causes a net mass flux 
along the slab axis. For stronger activity and/or 
larger system sizes, one can also find sponta- 
neously "shear-banded" flows in which succes- 
sive layers of material have very different shear 
rates. 

In a 2D (thin film) flow geometry, the pat- 
terns differed significantly from the ID case and 
were in good qualitative agreement with obser- 
vations of flow patterns in, for instance, con- 
centrated B. subtilus and Salmonella suspen- 
sions [191 EH [2H 122] ■ Fig. 2 shows some of the 
patterns which were found in 2D. In the apolar 
case, for values of £ just above the threshold, the 
instability first evolves into steady state "con- 
vective" rolls (Figs, la-b) [T5] . In some cases 
we also observed active bands or a succession 
of rolls (Fig. 2c), whereas deeper in the ac- 
tive phase we found that an initial array of rolls 
breaks up into what looks chaotic flow ( "bacte- 
rial turbulence" ) at low Reynolds number (Fig. 
2d-f). The equations of motion for polar active 
suspensions (see Appendix) can be treated with 
the same method. This is somewhat similar, 
in that above a threshold there is a transition 
to a spontaneously flowing state. In this case, 
however, we do not observe rolls or bands and 
the system jumps directly into the "turbulent" 
flowing state (see Fig. 2f-h, with 2h showing 
trajectories of tracer particles which highlight 
the chaotic appearance of the flow). 
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LB simulations have also proved helpful to char- 
acterise the rheology of active fluids. In Ref. [27] 
it was suggested that the viscosity of a con- 
tractile active fluid should diverge at the pas- 
sive isotropic-nematic transition in 2D. Simula- 
tions have confirmed and generalised this to the 
case of a 3D order parameter (albeit constrained 
to undergo a ID flow) where it was clarified 
that the divergence there occurs at the spinodal 
point (which unlike the 2D case is different from 
the isotropic-nematic transition point). LB sim- 
ulations also showed that upon increasing the 
density towards this spinodal point, one should 
observe a decrease in the viscosity for extensile 
fluids. 

The non-linear rheology should show even more 
striking behaviours [7J. Contractile fluids 
should strongly shear thicken for small forcing, 
and the extent of thickening should depend on 
the distance from the isotropic-nematic tran- 
sition. Within a bulk nematic phase, a for- 
mal yield stress (a stress threshold below which 
there is no flow in steady state) is predicted, 
whereas for larger forcing, these materials ex- 
hibit shear thinning and approach the unen- 
hanced (passive) viscosity at very large shear 
rates. 

In contrast, isotropic suspensions of concen- 
trated extensile particles (bacteria) should start 
from a low viscosity and thicken to again ap- 
proach the passive behaviour upon increasing 
the shear rate. Strikingly, a bulk oriented (ne- 
matic) phase of extensile active particles should 
show a zero effective viscosity for shear rates 
below some critical value. This is because the 
3D ordered system in the absence of stress (but 
constrained to have ID flow) spontaneously or- 
ganizes into two shear bands with flow in oppo- 
site directions, so that there is zero net velocity 
of the fluid. A finite relative velocity of the con- 
fining walls can now be accommodated, still at 
zero stress, by adjusting the relative amounts of 
the two bands. (Since there is now a finite shear 
rate and no shear stress, the viscosity is formally 
zero.) Finally, our studies suggested that when 
subjected to shear flow, extensile fluids should 
form shear bands more readily than passive ones 
close to the isotropic-nematic transition, while 
contractile fluids should shear-band less readily. 



These rheological predictions remain provi- 
sional, assuming as they do a ID flow profile 
and the results may well be modified when this 
assumption is relaxed. (Such work is now under- 
way in our group.) Nonetheless, the recent ad- 
vances in experimental techniques, which have 
made it possible to grow thin concentrated films 
of, e.g. B. subtilis, should render several of 
our results, such as the 2D predictions of Fig.2, 
testable in the near future. 

3 Blue phase simulations 

The "blue phases" (BPs) of chiral nemato- 
genic molecules offer spectacular examples of 
functional soft matter; each comprises a self- 
sustained network of disclinations [TU] embed- 
ded within a nematic matrix. (A disclination is 
a topological defect line, defined such that the 
nematic director rotates through a half-turn on 
treversing any circuit that encloses this line.) 
At high temperature, a fluid of chiral ne- 
matogenic molecules remains in the "isotropic 
phase", with no preferred orientation of the 
molecular axes. Upon cooling down the sam- 
ple, the molecules become oriented (Q is fi- 
nite), but due to molecular chirality the direc- 
tor field n rotates with spatial position, describ- 
ing everywhere a helix with a well-defined axis. 
(This is called the cholesteric phase.) However, 
very close to the transition, it is more advan- 
tageous locally for the director field to rotate 
in a helical fashion about any axis perpendic- 
ular to a straight line - this complicated pat- 
tern was named a "double twist cylinder" |10j . 
Mathematically, it is impossible to patch to- 
gether such double twist cylinders without cre- 
ating defects in between, and this frustration 
gives rise to the disclination network observed 
in the blue phases, and responsible for many 
of their remarkable physical properties. Most 
striking among these are their optical proper- 
ties (such phases can be made in all colours, 
not just blue) which stem from the presence of 
a lattice of disclinations with a unit cell whose 
size is comparable to the wavelength of visible 
light. In BP I and II (the two most common) 
this lattice has long range cubic order; BP III 
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however is certainly not cubic, and probably not 
ordered. 

BPs have a fascinating scientific history. They 
were first reported in the late 19th century, by a 
scientist named Reinitzer, and then long forgot- 
ten, until some new experimental interest arose 
in the 1960s and 1970s. Initial theories of BPs 
only came out in the 1980s (see e.g. [31]), when 
the concept of double twist was first proposed. 
At this stage BPs were widely considered to 
be of purely academic interest, mainly because 
they were only stable in a very narrow temper- 
ature range (about 1 K) close to the isotropic- 
cholesteric transition. In 1983, a world-leading 
expert on liquid crystals, F. C. Frank, said: 
"They [Blue phases] are totally useless, I think, 
except for one important intellectual use, that of 
providing tangible examples of topological odd- 
ities, and so helping to bring topology into the 
public domain of science, from being the pri- 
vate preserve of a few abstract mathematicians 
and particle theorists." [TO]. In the first decade 
of the 21st century, this view rather suddenly 
changed, following advances in fabrication that 
enormously increased the stability range of BPs, 
up to about 50 K [321 [33]. In May 2008 Sam- 
sung presented the first blue-phase based liquid 
crystal display at the annual SID International 
Symposium. This new display is able to oper- 
ate at a high frame frequency (240 Hertz), does 
not require costly alignment treatment at the 
boundaries of the liquid crystal, and may one 
day supersede current LC (twisted nematic) dis- 
play technologies. 

Traditional theories of BPs (see e.g. [34] ) 
were based on semi-analytical approximations. 
While they were extremely useful to gain a qual- 
itative understanding of the physics of BPs, 
these approaches had to rely on severe assump- 
tions for progress to be possible. For instance, 
when estimating the phase diagram, the Q ten- 
sor was approximated by a truncated Fourier 
series. Computational constraints only allowed 
very few harmonics to be considered. As a re- 
sult, when important experimental observations 
were mispredicted by the theory, it was not clear 
whether this was a drawback of the underlying 
free energy functional (the Landau-de Gennes 
free energy, see Appendix), or simply an arti- 



fact of the simplifications which were employed. 
These early theoretical works left unexplained 
the detailed shape of the phase diagram - BPI 
and BPII appeared in the wrong order on in- 
creasing the chirality k. (This parameter is 
proportional to the inverse helical pitch of the 
cholesteric; see Appendix.) With an electric 
field applied, the analytic theories were also un- 
able to account for the anomalous field-induced 
distortion (electrostriction) of BPI, nor to ex- 
plain why a new phase, named BPX, should be 
stable at all. Finally, most theories assume cu- 
bic symmetry, so cannot describe BPIII, or the 
"blue fog" , which is thought to comprise a net- 
work of disclination lines without long-range or- 
der. 

In recent years, LB simulations have been re- 
markably useful in filling most of these concep- 
tual theoretical gaps, and have significantly ex- 
tended our quantitative understanding of the 
physics of blue phases. Firstly, in Refs. [33 |3B] , 
Eq. [JJ was solved by non-hybrid LB in the 
absence of fluid flow (u = 0) with a set of 
initial conditions suggested from analytical ex- 
pressions for the infinite chirality limit of BPI, 
BPII and O5. (O5 is another disclination lattice, 
which was proposed as stable by early theories, 
but not observed in experiments.) This pro- 
cedure amounts to a free energy minimisation 
with a topological constraint; that is, the purely 
relaxational dynamics of Eq. [1] was shown to 
maintain the point group of the disclination net- 
work chosen initially. Therefore, this approach 
can be used to map out the full phase diagram as 
a function of chirality, k, and a reduced temper- 
ature parameter r (see Appendix for the math- 
ematical definition of these paramaters). This 
approach has the advantage of not making any 
approximations beyond those implicit in the se- 
lected (Landau-de Gennes) free energy. 
These LB simulations showed that the phase di- 
agram predicted by the continuum theory is ac- 
tually in good qualitative agreement with the 
experiments: BPI and BPII show up in the 
right order on increasing chirality, and O5 is 
relegated to unphysical regions in the phase 
diagram [35]. The structure of these phases 
is shown in Fig.3a,b where a surface is drawn 
around each defect line at a certain contour of 
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the ordering strength. (This creates a rendering 
of each disclination fattened tube.) The nu- 
merical phase diagram (Fig. 3c) is in good agree- 
ment with the experimental one (Fig. 3d). An 
interesting calculation in [35] has shown that 
even including an extra set of spherical har- 
monic in the analytical scheme of Ref. [31] is 
unable to reproduce the details shown by the 
numerics. This is an example in which simula- 
tions are extremely important to accurately find 
what the predictions of a given theory are. 
Similarly, LB simulations were performed in the 
presence of an electric field in Ref. [36]. It 
was found that under a small electric field the 
unit cell of BPII tended to elongate along the 
field direction, and shrink perpendicularly to it, 
whereas BPI displayed an opposite behaviour, 
once more in agreement with experiments. An 
intermediate field also turned the disclination 
network of BPI into a different structure (which 
it is tempting to identify with the experimen- 
tally observed BPX). Therefore, also when an 
electric field is present, the Landau-de Gennes 
free energy works remarkably well, although an- 
alytically tractable approximations of the re- 
sulting equations are not adequate to capture 
its predictions. 

As well as allowing one to find equilibrium states 
under different fields and thermodynamic con- 
ditions, LB of course comes into its own for dy- 
namical problems in which fluid flow cannot be 
ignored. For instance, the existence of a discli- 
nation network affects the response of a BP to 
an imposed Poiseuille flow [39J. Here, LB work 
has shown that flow can lead to the unzipping of 
disclinations of integer topological charge. (On 
a circuit around such a defect line, the direc- 
tor rotates through a whole number of turns, 
rather than the half-turn around a standard 
disclination; these higher-order lines can form 
metastable networks in the absence of an ap- 
plied flow.) Flow also can cause the bending 
and twisting of the BPI and BPII disclination 
networks. This bending and twisting lead to an 
elastic component in the rheological behaviour, 
and result the simulations predict "perme- 
ative" flows, in which the molecules comprising 
the blue phase flow through a static disclination 
pattern whose geometry is hardly perturbed by 



the underlying molecular transport. The same 
kind of flow also occurs when cholesterics are 
sheared by small forcing along the direction of 
their helical axis [ID]. Typically, BPs also dis- 
play significantly shear thinning behaviour, as 
a strong enough flow disrupts the disclination 
network in a manner that reduces the elastic 
stress. 

Until very recently, LB work on BPs was lim- 
ited to one unit cell of the disclination lat- 
tice, within which several disclination cores are 
present and require a fine enough discretisation 
to be correctly resolved. However, supercom- 
puters now allow supra-unit cell simulations of 
BPs, the first account of which we have recently 
given in Ref. [H], where we studied the do- 
main growth dynamics of a BPII domain inside 
a cholesteric or isotropic "slab" , in a parameter 
region in which the BP is the thermodynam- 
ically stable state. The simulations give evi- 
dence of an intriguing domain growth kinetics. 
For small values of the chirality, the growth is 
slow and the resulting blue phase has no or few 
defects. When the chirality exceeds a certain 
threshold, however, the advancing disclination 
network changes its symmetry and reconstruct 
into a new hexagonal phase, which is so far un- 
documented in experiments. (This process is 
shown in Fig. 3e,f.) It would be of interest to 
determine whether this new BP is a metastable 
structure found due to the geometry we have fo- 
cussed on (we considered just one planar slice of 
unit cells) or has a wider physical meaning. In 
all cases, these simulations are encouraging as 
they suggest that large scale simulations of BPs 
are within reach. This is of course of interest 
to the modelling of real devices, which can be 
manufactured in the micron scale which we can 
consider computationally. 

4 Conclusions and future 
prospects 

We hope that this selection of results has shown 
that LB simulations of liquid crystals are poten- 
tially extremely powerful in gaining new insights 
into the physics which is contained in the hy- 
drodynamic equations of motion of liquid crys- 
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talline fluids. For both active fluids and blue 
phases, it would have been very difficult to com- 
pare theory and experiments - even at a qual- 
itative level - without using these simulations. 
Although such comparisons remain in their in- 
fancy for active nematics, the spontaneous flow 
patterns in (for instance) concentrated bacte- 
rial suspensions, simulated via the continuum 
equations of motion of an active liquid crystal, 
are qualitatively comparable with the patterns 
observed in the experiments. In the case of 
blue phases, such comparisons are more clear 
cut. Here it was thought that the classic theory 
based on a Landau-de Gennes free energy was 
missing some physics because (for example) the 
phase diagram was poorly predicted and anoma- 
lous electrostriction in BPI was not found. Re- 
markably, LB has shown that this was a draw- 
back of the approximations used to make ana- 
lytical progress, and not of the original theory, 
which is qualitatively and semi-quantitatively 
accurate. 

The fields we have covered in this short review 
are, of course, still full of open questions, and 
we hope that future LB simulations will play an 
important role in clarifying some of these. 
In active fluids, it will be important to 
characterise the flow patterns in fully three- 
dimensional active suspensions, and also to ex- 
tend the treatment we have covered here to the 
case in which there are density fluctuations or 
inhomogeneity in the fluid. Another related is- 
sue would be to study "active emulsions", in 
which droplets of active gels are suspended in 
an aqueous passive medium, possibly enclosed 
by an elastic membrane. Ultimately, it would 
be very exciting if continuum theories like the 
one we solved numerically may be applied to, 
for instance, suspensions of cell extracts in an 
extracellular matrix. From the theory point of 
view, it appears that an urgent issue is to clar- 
ify to what extent active fluids faithfully rep- 
resent concentrated suspensions of motile parti- 
cles, or swimmers, by, for instance, comparing 
the results of continuum simulations to those 
of more microscopic models with fully resolved 
swimmers, which can also be treated via LB 
(though of a different kind than the one pre- 
sented here) [221 S21 S3] - 



Large scale simulations of blue phases will also 
be likely to be important in the future. From 
an application point of view, the exciting po- 
tential of BP devices can ultimately be fully 
exploited if we manage to reach a quantitative 
understanding of their thermodynamics, their 
switching dynamics, and the role of flow. Supra- 
unit cell simulations are needed to this end, be- 
cause the field leads to unit cell deformations 
and may cause full scale reconstruction of the 
disclination network. From a more fundamen- 
tal point of view, we do not have a satisfactory 
understanding of non-cubic blue phases. Most 
notably, the structure of BPIII - the "blue fog" 
- is still not understood to date, and we hope 
that large scale simulations of amorphous discli- 
nation networks may shed some light on this 
elusive problem. 

We are grateful to G. P. Alexander, S. M. 
Fielding, A. N. Morozov, E. Orlandini and J. 
M. Yeomans for useful discussions. We ac- 
knowledge EPSRC grants EP/E045316/1 and 
EP/E030 173/1 for funding, and computer time 
on Hector funded by EP/F054750/1. MEC 
holds a Royal Society Research Professorship. 
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Fig. 1 (a) Turbulence in a sessile droplet of 
B. subtilis, viewed from below a petri 
dish. The horizontal line is the edge 
of the droplet (picture taken from 
Fig. 3 of Ref. |2Q]). The scale bar 
is 35 /jm. (b) Flow pattern in a sim- 
ilar bacterial droplet - the arrow at 
the right stands for a speed of 35 
//m/s (picture taken from Fig. 4 of 
Ref. (2D]). 

Fig. 2 Selected results from active fluid 
simulations. We only plot the ve- 
locity field, resulting from LB so- 
lutions of Eq. [2j The top row (a- 
c) shows stationary states obtained 
for apolar extensile fluids with mod- 
erate activity: it can be seen that 
the spontaneous flow has the shape 
of rolls (a,b) or of bands, in general 
tilted (c). The middle row show non- 
stationary "turbulent" solutions for 
larger values of the activity. The 
bottom row shows solution of the 
equations of motion of polar active 
gels. In (g) there is no self-advection 
term, so that the fluid is equivalent 
to an apolar gel, whereas in (h-i) 
this term is switched on. In (i) we 
plot the trajectories of 3 tracer par- 
ticles, which show the "turbulent" 
nature of the flow. Parameters com- 
mon to the apolar runs in (a)-(f) are: 
7 = 3 (ensuring that we work in the 
ordered phase), £ = 0.7 (which to- 
gether with our choice of 7 selects 
flow-aligning liquid crystals), T = 
0.33775, K = 0.08, r) = 0.57. The 
activity parameter ( was 0.001 (a), 
0.002 (b), 0.01 (c), and 0.04 (d-f). 
Parameters for the polar runs in (g- 
i) are ( = 0.02, A = 1.1, K = 0.04, 
T = 0.3, 77 = 1.67, and w = (g), or 
w = 0.01 (h,i). 

Fig. 3 The top row shows the disclination 
lattices of BPI (a) and BPII (b), 
as obtained from LB simulations - 
2 unit cells in each directions are 
shown. The second row shows a 



computational (c, within the one 
elastic constant approximation), and 
an experimental (d) typical phase di- 
agram for blue phases (picture taken 
from Fig. 3 of Ref. [H7]). The key 
feature is that LB simulations pre- 
dict the correct order of appearance 
of BPI and BPII upon increasing 
the chirality. Note that BPIII is 
not a cubic phase hence is not in- 
cluded in the theoretical phase di- 
agram. The bottom two rows (e- 
f) show dynamical states obtained 
when a BPII domain grows inside 
an initially cholesteric matrix (see 
Ref. [2]). The reduced tempera- 
ture was r = and the chirality was 
k = 2 (e) or k = 1 (f). 
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Figure 1: 
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Figure 3: 
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Appendix: Hydrodynamic equations 
of motion for active and passive liquid 
crystalline fluids 

In this Appendix we review the equations of mo- 
tion for (active and passive) liquid crystalline 
fluids, which we solve by lattice Boltzmann sim- 
ulations. These are the equations used to gen- 
erate the results reviewed in our work. 
We first describe the thermodynamics of a liquid 
crystalline fluid in the absence of active stresses. 
This covers cholesterics and blue phases (and 
also active gels within a passive phase). We em- 
ploy a Landau-de Gennes free energy J 7 , whose 
density we indicate by /. The free energy den- 
sity can be written as a sum of two contribu- 
tions, fx and /2- The first is a bulk contribu- 
tion, 



+ 



Ao ( . 
2 1 



)Ql P 



A)7 



while the second is a distortion term. For 
nonchiral liquid crystals, we take the (standard) 
one elastic constant approximation [9] 



h 



Where Aq is a constant, 7 controls the magni- 
tude of order (it may be viewed as an effective 
temperature or concentration for thermotropic 
and lyotropic liquid crystals respectively), while 
K is an elastic constant. To describe cholester- 
ics, we employ the slightly generalised distortion 
free energy, which is again standard [TP] : 

K r 



a/3, 



Here and in what follows Greek indices denote 
cartesian components and summation over re- 
peated indices is implied. 

For blue phases, it is customary to identify the 
position in thermodynamic parameter space via 
the chirality, k, and the reduced temperature, 
r. These may be defined in terms of previous 
quantities via 



108Kqi 
\l 4>7 
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' 1-7/3 ' 
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(Note that the reduced temperature was defined 



in older literature as 



T 
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When needed, the anchoring of the director field 
on the boundary surfaces (Fig. 2) to a chosen 
unit vector n° is ensured by adding a surface 
term in the free energy density 



fs 

Qa/3 



2^o(Q a /3 



S (n a n^ - 6 a p/3) 



The parameter Wq controls the strength of the 
anchoring, while So determines the degree of the 
surface order. If the surface order is to equal the 
bulk order, S should be set equal to q, the or- 
der parameter in the bulk (3/2 times the largest 
eigenvalue of the Q tensor) . Wo is large (strong 
anchoring) in what follows. 
The equation of motion for Q is taken to be [9J 

(dt + u- V)Q - S(W, Q) = TH + AQ 

where T is a collective rotational diffusion con- 
stant, and A is an activity parameter which for 
simplicity we set to zero in our simulations. 
(The resulting term can anyway be absorbed 
into a shift of A and/or 7 in the free energy.) 
The first term on the left-hand side of the equa- 
tion above is the material derivative describing 
the usual time dependence of a scalar quantity 
advected by a fluid with velocity u. This is mod- 
ified for rod-like molecules by a second term 

S(W,Q) = (£D + w)(Q + I/3) 
+ (Q + I/3)(£D-w) 
- 2£(Q + I/3)Tr(QW) 

where Tr denotes the tensorial trace, while 
D = (W + W T )/2 and u = (W - W T )/2 
are the symmetric part and the anti-symmetric 
part respectively of the velocity gradient tensor 
W a p = dpu a . The constant £ depends on the 
molecular details of a given liquid crystal, and 
determines, together with 7, whether a liquid 
crystal is flow aligning or flow tumbling (we re- 
strict to the former case in this work) . The first 
term on the right-hand side of the order param- 
eter evolution equation describes the relaxation 
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of the order parameter towards the minimum of 
the free energy. The molecular field H which 
provides the driving motion is given by 

The fluid velocity, u, obeys the continuity equa- 
tion and the Navier-Stokes equation, whose in- 
compressible limit is Eq. [2] in which Ti a p = 
■Qpassive _|_ reactive _ rj ne s ^ Tess tensor U^ slve nec- 

essary to describe ordinary LC hydrodynamics 
is (up to an isotropic pressure term) given by: 



T-rpassive 
U a,fl 



5F 



whereas the active term is given by, in leading 
order 

where ( is an activity constant [17]. Note that 
with the sign convention chosen here ( > cor- 
responds to extensile rods and ( < to contrac- 
tile ones [TTj . 

We can also use a variant of these equations to 
study polar active gels. The order parameter is 
this time a vector P a (with variable magnitude) 
as there is no longer head-tail symmetry in the 
system. The equations of motion we used in 
our LB simulations (reported in Fig. 2, bottom 
row) , are a simplified version of those presented 
in [H]. The equation governing the evolution of 
the vectorial order parameter is 



The "molecular field" is now given by h a = 
—8J r P oi/5p a where J-p i is the free energy for a 
polar active nematic, whose density is (see also 
Ref. [H] where a more general form is used): 



/ 



4 



P\* + K(d a P p 



(3) 



where K is an elastic constant, b > 0, and we 
chose a = —b to ensure that the minimum of 
the free energy is with |P| = 1. (Note that in 
this model P this time denotes a vector rather 
than a tensor). 

The Navier-Stokes equation is as in the tensorial 
model, but the stress tensor this time is 



2 (Pah 



PrK 



\ (Pah 



(P a Ps 



As mentioned in the text, the active term in the 
stress tensor, proportional to (, has therefore 
the same form for apolar and polar active gels. 
A few limits of the two theories considered 
above are worth noting. The tensorial model 
we have written down is equal to, for ( = 0, the 
Beris-Edwards model for liquid crystal hydrody- 
namics. Analogously, for ( = w = the polar 
model reduces to the Leslie-Ericksen model of 
nematodynamics. For w — 0, and a sample of 
uniaxial active liquid crystals, with a spatially 
uniform degree of orientational order, the ten- 
sorial model may be mapped onto the vectorial 
one (see [13] for a proof of this). 



[d t + (up + wPp) dp)} P c 



In this equation, w is another active term, due 
to swimming, which causes self-advection of the 
order parameter, while A is a material depen- 
dent constant - positive for rod-like molecules. 
If | A | > 1 the liquid crystalline passive phase 
is flow- aligning, otherwise it is flow-tumbling. 
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